#ifndef ppform_h
#define ppform_h

#include<iostream>
#include<cmath>
#include<vector>
#include<cstring>
#include<fstream>
#include"Eigen/Dense"

using namespace std;
using namespace Eigen;

const double _epsL = 10 * std::numeric_limits<double>::epsilon();
int flag=1;
class Function{
public:
    virtual double operator () (const double &x) const = 0;
    virtual double diff (const double &x) const{
        return ((*this)(x+_epsL)-(*this)(x-_epsL)) / (2*_epsL);
    }
    virtual double diff2(const double &x) const{
        return ((*this)(x+2*_epsL)+(*this)(x-2*_epsL)-2*(*this)(x))/(4*_epsL);
    }
};

/*class EquationSolver{
public:
    virtual double solve() = 0;
};
*/
class Point{
public:
    vector<double> x;//index from 1, means x_0 is the first elements
    vector<double> f;
    Point(){};
    Point(vector<double> x,vector<double> f):f(f),x(x){};
    bool equality(Point X){   //determine whether two or more same x's map to different y's;
        for(int i=0;i<X.x.size();i++){
            for (int j=0;j<X.f.size();j++){
                if(x[i]==x[j]&&f[i]!=f[j]){
                    return false;
                    break;
                }
            }
        }
        return true;
    }
    
    bool determine(Point X){
        if(X.x.size()!=X.f.size()){
            cout<<"ERROR";
            return false;
        }
        else if(equality(X)==false){
            return false;
        }
        else return true;
    }
    void showpoint(Point X){
        if(determine(X)==true&&equality(X)==true){
            for(int i=0; i<X.x.size(); i++)
            {
                cout<<X.x[i]<<" "<<X.f[i]<<endl;
            }
        }
        else
            cout<<"ERROR";
    }
    double getsize(Point X){
        if(determine(X)==true&&equality(X)==true){
            return X.x.size();}
        else{
            return -1;
        }
    }
};

class PiecewisePoly:public Point{
private:
    Function &Func;
    Point X;
    int n=X.getsize(X);
    vector<double> lamda,miu,h,d,ff,diff,delta;
    vector<double> mp(){
        vector<double> m;
        for(int i=0;i<X.x.size();i++){
            m.push_back((X.x[i]+X.x[i+1])/2.0);
        }
        return m;
    }
   // double ANS[6000][4];
    void find_h(Point X){
        for (int i=0;i<n;i++){
            h.push_back(X.x[i+1]-X.x[i]); 
        }
    }
    MatrixXf solve_axb(MatrixXf A,MatrixXf b){
        MatrixXf x;
        x=A.fullPivLu().solve(b);
        return x;
    }
    void lamda_miu(vector<double> x){
        for(int i=1;i<size(x);i++){
            lamda.push_back((x[i+1]-x[i])/(x[i+1]-x[i-1]));
            miu.push_back((x[i]-x[i-1])/(x[i+1]-x[i-1]));
        }
    }
   
public:
    PiecewisePoly(Function &Func,Point X):Func(Func),X(X){};
    /*void diff_h_delta(Point X,double xx){
        for(int i=0;i<X.x.size();i++){
            diff.push_back((X.f[i+1]-X.f[i])/(X.x[i+1]-X.x[i]));
            h.push_back(X.x[i+1]-X.x[i]);
            delta.push_back(xx-X.x[i]);
    
        }
    }*/
    void condition(string z,double xx){
        if(z=="complete"){
            MatrixXf A(n,n);
            MatrixXf b(n,1);
            MatrixXf M(X.x.size(),1);
            for(int i=0;i<n;i++){
                if(i!=0&&i!=n-1){
                    b(i,0)=6*((X.f[i+1]-X.f[i])/(X.x[i+1]-X.x[i])-(X.f[i]-X.f[i-1])/(X.x[i]-X.x[i-1]))/(X.x[i+1]-X.x[i-1]);
                }
                b(0,0)=6*(((X.f[1]-X.f[0])/(X.x[1]-X.x[0]))-Func.diff(X.x[0]))/(X.x[1]-X.x[0]);
                b(n-1,0)=6*((-(X.f[n-1]-X.f[n-2])/(X.x[n-1]-X.x[n-2]))+Func.diff(X.x[n-1]))/(X.x[n-1]-X.x[n-2]);
            }
            for(int i=0;i<n;i++){
                for(int j=0;j<n;j++){
                    if(i==j){
                        A(i,j)=2;
                    }
                    else if(i>0&&i<n-1){
                        A(i,i+1)=(X.x[i+1]-X.x[i])/(X.x[i+1]-X.x[i-1]);
                        A(i,i-1)=(X.x[i]-X.x[i-1])/(X.x[i+1]-X.x[i-1]);
                    }
                    A(0,1)=1;
                    A(n-1,n-2)=1;
                }
            }
            M=solve_axb(A, b);
            /*cout<<A<<endl;
            cout<<b<<endl;
            cout<<M<<endl;*/
            Point Z;
            int count=0;
            for(int i=0;i<X.x.size()-1;i++){
                if(xx>=X.x[i]&&xx<=X.x[i+1]){
                    Z.x.push_back(xx);
                    Z.f.push_back(X.f[i]+((X.f[i+1]-X.f[i])/(X.x[i+1]-X.x[i])-1.0*(2*M(i,0)+M(i+1,0))*(X.x[i+1]-X.x[i])/6.0)*(xx-X.x[i])+M(i,0)*pow(xx-X.x[i],2)/2.0+(M(i+1,0)-M(i,0))*pow(xx-X.x[i],3)/(6.0*(X.x[i+1]-X.x[i])));
                    ofstream out_to_file("A.txt");
                    out_to_file<<Z.x[count]<<" "<<Z.f[count]<<endl;
                    cout<<Z.x[count]<<" "<<Z.f[count]<<endl;
                    count++;
                    out_to_file.close();
                }
                else continue;
            }
            
            Point Error;
            vector<double> midpoint;
            midpoint=mp();
            int cnt=0;
            //cout<<"Below are plenty of points to approximate the given functions:"<<endl;
            if (flag!=0){
                for(int i=0;i<X.x.size()-1;i++){
                    if(midpoint[cnt]>=X.x[i]&&midpoint[cnt]<=X.x[i+1]){
                        Error.x.push_back(midpoint[cnt]);
                        Error.f.push_back(abs(Func(midpoint[cnt])-(X.f[i]+((X.f[i+1]-X.f[i])/(X.x[i+1]-X.x[i])-1.0*(2*M(i,0)+M(i+1,0))*(X.x[i+1]-X.x[i])/6.0)*(midpoint[cnt]-X.x[i])+M(i,0)*pow(midpoint[cnt]-X.x[i],2)/2.0+(M(i+1,0)-M(i,0))*pow(midpoint[cnt]-X.x[i],3)/(6.0*(X.x[i+1]-X.x[i])))));
                        
                        ofstream out_to_file("A_Error.txt");
                        out_to_file<<Error.x[cnt]<<" "<<Error.f[cnt]<<endl;
                        //cout<<Error.x[cnt]<<" "<<Error.f[cnt]<<endl;
        
                        cnt++;
                        //out_to_file.close();
                    }
                    else continue;
                }
            }
            flag=0;
           /* for(int i=0;i<cnt;i++){
            
                cout<<Error.x[i]<<" "<<Error.f[i]<<endl;
            }*/
        }
        else if(z=="natural"){
            MatrixXf A(n,n);
            MatrixXf b(n,1);
            MatrixXf M(n,1);
            for(int i=0;i<n;i++){
                if(i!=0&&i!=n-1){
                    b(i,0)=6*((X.f[i+1]-X.f[i])/(X.x[i+1]-X.x[i])-(X.f[i]-X.f[i-1])/(X.x[i]-X.x[i-1]))/(X.x[i+1]-X.x[i-1]);
                }
                else b(i,0)=0;
            }
            for(int i=0;i<n;i++){
                for(int j=0;j<n;j++){
                    if((i==0&&j==0)||(i==n-1&&j==n-1)){
                        A(i,j)=1;
                    }
                    else if(i==j){
                        A(i,j)=2;
                    }
                    else if(i>0&&i<n-1){
                        A(i,i-1)=(X.x[i+1]-X.x[i])/(X.x[i+1]-X.x[i-1]);
                        A(i,i+1)=(X.x[i]-X.x[i-1])/(X.x[i+1]-X.x[i-1]);
                    }
                }
            }
            M=solve_axb(A,b);
            /*cout<<A<<endl;
            cout<<b<<endl;
            cout<<M;*/
            Point Z;
            int count=0;
            for(int i=0;i<X.x.size()-1;i++){
                if(xx>=X.x[i]&&xx<=X.x[i+1]){
                    Z.x.push_back(xx);
                    Z.f.push_back(X.f[i]+((X.f[i+1]-X.f[i])/(X.x[i+1]-X.x[i])-1.0*(2*M(i,0)+M(i+1,0))*(X.x[i+1]-X.x[i])/6.0)*(xx-X.x[i])+M(i,0)*pow(xx-X.x[i],2)/2.0+(M(i+1,0)-M(i,0))*pow(xx-X.x[i],3)/(6.0*(X.x[i+1]-X.x[i])));
                    ofstream out_to_file("A.txt");
                    out_to_file<<Z.x[count]<<" "<<Z.f[count]<<endl;
                    cout<<Z.x[count]<<" "<<Z.f[count]<<endl;
                    count++;
                    out_to_file.close();
                }
                else continue;
            }
            Point Error;
            vector<double> midpoint;
            midpoint=mp();
            int cnt=0;
            //cout<<"Below are plenty of points to approximate the given functions:"<<endl;
            if (flag!=0){
                for(int i=0;i<X.x.size()-1;i++){
                    if(midpoint[cnt]>=X.x[i]&&midpoint[cnt]<=X.x[i+1]){
                        Error.x.push_back(midpoint[cnt]);
                        Error.f.push_back(abs(Func(midpoint[cnt])-(X.f[i]+((X.f[i+1]-X.f[i])/(X.x[i+1]-X.x[i])-1.0*(2*M(i,0)+M(i+1,0))*(X.x[i+1]-X.x[i])/6.0)*(midpoint[cnt]-X.x[i])+M(i,0)*pow(midpoint[cnt]-X.x[i],2)/2.0+(M(i+1,0)-M(i,0))*pow(midpoint[cnt]-X.x[i],3)/(6.0*(X.x[i+1]-X.x[i])))));
                        
                        ofstream out_to_file("A_Error.txt");
                        out_to_file<<Error.x[cnt]<<" "<<Error.f[cnt]<<endl;
                       // cout<<Error.x[cnt]<<" "<<Error.f[cnt]<<endl;
        
                        cnt++;
                        //out_to_file.close();
                    }
                    else continue;
                }
            }
            flag=0;
           /* for(int i=0;i<cnt;i++){
            
                cout<<Error.x[i]<<" "<<Error.f[i]<<endl;
            }*/
        }
        else if(z=="notaknot"){
            MatrixXf A(n,n);
            MatrixXf b(n,1);
            MatrixXf M(n,1);
            for(int i=0;i<n;i++){
                if(i!=0&&i!=n-1){
                    b(i,0)=6*((X.f[i+1]-X.f[i])/(X.x[i+1]-X.x[i])-(X.f[i]-X.f[i-1])/(X.x[i]-X.x[i-1]))/(X.x[i+1]-X.x[i-1]);
                }
                else b(i,0)=0;
            }
            for(int i=0;i<n;i++){
                for(int j=0;j<n;j++){
                    if(i==0){
                        A(i,0)=-(X.x[2]-X.x[1])/(X.x[2]-X.x[0]);
                        A(i,1)=1;
                        A(i,2)=-(X.x[1]-X.x[0])/(X.x[2]-X.x[0]);
                    }
                    else if(i==n-1){
                        A(i,n-3)=-(X.x[i]-X.x[i-1])/(X.x[i]-X.x[i-2]);
                        A(i,n-2)=1;
                        A(i,n-1)=-(X.x[i-1]-X.x[i-2])/(X.x[i]-X.x[i-2]);
                    }
                    
                    else if(i==j){
                        A(i,j)=2;
                    }
                    else if(i>0&&i<n-1){
                        A(i,i-1)=(X.x[i+1]-X.x[i])/(X.x[i+1]-X.x[i-1]);
                        A(i,i+1)=(X.x[i]-X.x[i-1])/(X.x[i+1]-X.x[i-1]);
                    }
                }
            }
            M=solve_axb(A,b);
            /*cout<<A<<endl;
            cout<<b<<endl;
            cout<<M<<endl;*/
            
            Point Z;
            int count=0;
            for(int i=0;i<X.x.size()-1;i++){
                if(xx>=X.x[i]&&xx<=X.x[i+1]){
                    Z.x.push_back(xx);
                    Z.f.push_back(X.f[i]+((X.f[i+1]-X.f[i])/(X.x[i+1]-X.x[i])-1.0*(2*M(i,0)+M(i+1,0))*(X.x[i+1]-X.x[i])/6.0)*(xx-X.x[i])+M(i,0)*pow(xx-X.x[i],2)/2.0+(M(i+1,0)-M(i,0))*pow(xx-X.x[i],3)/(6.0*(X.x[i+1]-X.x[i])));
                    ofstream out_to_file("A.txt");
                    out_to_file<<Z.x[count]<<" "<<Z.f[count]<<endl;
                    cout<<Z.x[count]<<" "<<Z.f[count]<<endl;
                    count++;
                    out_to_file.close();
                }
                else continue;
            }
            Point Error;
            vector<double> midpoint;
            midpoint=mp();
            int cnt=0;
            //cout<<"Below are plenty of points to approximate the given functions:"<<endl;
            if (flag!=0){
                for(int i=0;i<X.x.size()-1;i++){
                    if(midpoint[cnt]>=X.x[i]&&midpoint[cnt]<=X.x[i+1]){
                        Error.x.push_back(midpoint[cnt]);
                        Error.f.push_back(abs(Func(midpoint[cnt])-(X.f[i]+((X.f[i+1]-X.f[i])/(X.x[i+1]-X.x[i])-1.0*(2*M(i,0)+M(i+1,0))*(X.x[i+1]-X.x[i])/6.0)*(midpoint[cnt]-X.x[i])+M(i,0)*pow(midpoint[cnt]-X.x[i],2)/2.0+(M(i+1,0)-M(i,0))*pow(midpoint[cnt]-X.x[i],3)/(6.0*(X.x[i+1]-X.x[i])))));
                        
                        ofstream out_to_file("A_Error.txt");
                        out_to_file<<Error.x[cnt]<<" "<<Error.f[cnt]<<endl;
                       // cout<<Error.x[cnt]<<" "<<Error.f[cnt]<<endl;
        
                        cnt++;
                        //out_to_file.close();
                    }
                    else continue;
                }
            }
            flag=0;
           /* for(int i=0;i<cnt;i++){
            
                cout<<Error.x[i]<<" "<<Error.f[i]<<endl;
            }*/
        }
        else if(z=="linear"){
            Point Z;
            int count=0;
            for(int i=0;i<X.x.size()-1;i++){
                if(xx>=X.x[i]&&xx<=X.x[i+1]){
                    Z.x.push_back(xx);
                    Z.f.push_back(((xx-X.x[i+1])/(X.x[i]-X.x[i+1]))*X.f[i]+((xx-X.x[i])/(X.x[i+1]-X.x[i]))*X.f[i+1]);
                    ofstream out_to_file("A.txt");
                    out_to_file<<Z.x[count]<<" "<<Z.f[count]<<endl;
                    cout<<Z.x[count]<<" "<<Z.f[count]<<endl;
                    count++;
                    out_to_file.close();
                }
                else continue;
            }
        }
        else cout<<"ERROR"<<endl;
    }
   
};

class Bspline:public Point{
private:
    Function &Func;
    Point X;
    int order;
    double x_0,x_m1,x_m2,x_Np1,x_Np2,x_Np3;
    int n=X.getsize(X);
    void get_newpoint(Point X){
        x_0=X.x[0]-(X.x[1]-X.x[0]);
        x_m1=x_0-(X.x[1]-X.x[0]);
        x_m2=x_m1-(X.x[1]-X.x[0]);
        x_Np1=X.x[n-1]+(X.x[1]-X.x[0]);
        x_Np2=x_Np1+(X.x[1]-X.x[0]);
        x_Np3=x_Np2+(X.x[1]-X.x[0]);
    }
    MatrixXf solve_axb(MatrixXf A,MatrixXf b){
        MatrixXf x;
        x=A.fullPivLu().solve(b);
        return x;
    }
public:
    Bspline(Function &Func,Point X,int order):Func(Func),X(X),order(order){};
    double B_i(double xx,int i,int order){
       if(order==0){
            if(xx>i-1&&xx<=i){
                return 1.0;
            }
            else{
                return 0;
            }
        }
        else{
            return (xx-i+1)/order*B_i(xx,i,order-1)+(i+order-xx)/order*B_i(xx,i+1,order-1);
        }
    };
    double B_n(int n, int i, double xx){
            if ( n==0){
                if (xx > X.x[i+1] && xx<= X.x[i+2]){
                    return 1.0;
                }
                else{
                    return 0;
                }
            }
            else{
            return (xx - X.x[i+1])/(X.x[i+n+1] - X.x[i+1])*B_n(n-1,i,xx) + (X.x[i+n+2] - xx)/(X.x[i+n+2] - X.x[i+2])*B_n( n-1,i+1,xx);
            }
        }
    void condition(string z,double xx,int i){
        if(z=="complete"){
            MatrixXf b(n+2,1);
            MatrixXf A(n+2,n+2);
            MatrixXf a(n+2,1);
            for(int j=1;j<n+1;j++){
                b(j,0)=X.f[j];
            }
            b(0,0)=Func.diff(X.x[0]);
            b(n+1,0)=Func.diff(X.x[n-1]);
            for(int j=1;j<n+1;j++){
                A(j,j-1)=B_n(X.x[j-1],j-2,3);
                A(j,j)=B_n(X.x[j-1],j-1,3);
                A(j,j+1)=B_n(X.x[j-1],j,3);
            }
            A(0,0)=-3*B_n(X.x[0],0,2)/(X.x[1]-(X.x[0]-2*(X.x[1]-X.x[0])));
            A(0,1)=3*B_n(X.x[0],0,2)/(X.x[1]-(X.x[0]-2*(X.x[1]-X.x[0])))-3*B_n(X.x[0],1,2)/(X.x[2]-(X.x[0]-(X.x[1]-X.x[0])));
            A(0,2)=3*B_n(X.x[0],1,2)/(X.x[1]-(X.x[0]-(X.x[1]-X.x[0])));
            A(n+1,n-1)=-3*B_n(X.x[n-1],n-1,2)/(X.x[n-1]+(X.x[1]-X.x[0])-X.x[n-3]);
            A(n+1,n)=3*B_n(X.x[n-1],n-1,2)/(X.x[n-1]+(X.x[1]-X.x[0])-X.x[n-3])-3*B_n(X.x[n-1],n,2)/(X.x[n-1]+2*(X.x[1]-X.x[0])-X.x[n-2]);
            a=solve_axb(A,b);
            cout<<A<<endl;
            cout<<b<<endl;
            cout<<a<<endl;
        }
        else if(z=="natural"){
            MatrixXf b(n+2,1);
            MatrixXf A(n+2,n+2);
            MatrixXf a(n+2,1);
            for(int i=0;i<n;i++){
                b(i,0)=X.f[i];
            }
            b(n,0)=0;
            b(n+1,0)=0;
        }
        else if(z=="notaknot"){
            MatrixXf b(n+2,1);
            MatrixXf A(n+2,n+2);
            MatrixXf a(n+2,1);
            for(int i=0;i<n;i++){
                b(i,0)=X.f[i];
            }
            b(n,0)=0;
            b(n+1,0)=0;
        }
        else if(z=="cardinal_quadratic"){
            Point Z;
            /*for(int k=0;k<X.x.size();k++){
                Z.x.push_back(xx);
                if(xx>i-1&&xx<=i){
                    Z.f.push_back((pow(xx-i+1,2))/2.0);
                }
                else if(xx>i&&xx<=i+1){
                    Z.f.push_back(0.75-pow(xx-i-0.5,2));
                }
                else if(xx>i+1&&xx<=i+2){
                    Z.f.push_back(pow(i+2-xx,2)/2.0);
                }
                else Z.f.push_back(0);
            }*/
        }
        else if(z=="cardinal_cubic"){
            Point Z;
            /*for(int k=0;k<X.x.size();k++){
                Z.x.push_back(xx);
                if(xx>i-1&&xx<=i){
                    Z.f.push_back((pow(xx-i+1,3))/6.0);
                }
                else if(xx>i&&xx<=i+1){
                    Z.f.push_back(2.0/3.0-pow(i+1-xx,2)*(xx-i+1)*0.5);
                }
                else if(xx>i+1&&xx<i+3){
                    Bspline(Func,X,order).condition("cardinal_cubic",2*i+2-xx,i);
                }
                else Z.f.push_back(0);
            }*/
        }
        else if(z=="linear"){
            Point Z;
            MatrixXf a(n,1);
            double result=0.0;
            for(int k=0;k<X.x.size();k++){
                a(k,0)=X.f[k];
                result=result+a(k,0)*B_n( xx, k, order);
            }
            Z.x.push_back(xx);
            Z.f.push_back(result);
        }
        else {
            cout<<"ERROR"<<endl;
            exit(-1);
        }
    }
    /*double getBspline_value(double x,int i,int k,vector<double> t){
        if(i<0){
            cout<<"ERROR"<<endl;
            return 0;
        }
        else{
            double value=0.0,value1=0.0,value2=0.0;
            if(k==0){
                if(x<t[i]||x>t[i+1]){
                    return value;
                }
                else{
                    value=1.0;
                    return value;
                }
            }
            else if(k>0){
                if(x<t[i]||x>t[i+k+1]){
                    return value;
                }
                else{
                    double a1=0.0,a2=0.0,deno=0.0;
                    deno=t[i+k]-t[i];
                    if(deno==0){
                        a1=0;
                    }
                    else{
                        a1=(x-t[i])/deno;
                    }
                    deno=t[i+k+1]-t[i+1];
                    if(deno==0){
                        a2=0;
                    }
                    else{
                        a2=(t[i+k+1]-x)/deno;
                        
                    }
                    value1=a1*getBspline_value(x, i, k-1, t);
                    value2=a2*getBspline_value(x, i+1, k-1, t);
                    value=value1+value2;
                }
                return value;
            }
            else{
                cout<<"ERROR"<<endl;
                return 0;
            }
        }
    }
    double getBspline_value_i(double x,int i,int k,vector<double> t){
        if(i<-1){
            cout<<"ERROR"<<endl;
            return 0;
        }
        else{
            double value=0.0,value1=0.0,value2=0.0;
            if(k==0){
                if(x<t[i+1]||x>t[i+2]){
                    return value;
                }
                else{
                    value=1.0;
                    return value;
                }
            }
            else if(k>0){
                if(x<t[i+1]||x>t[i+k+2]){
                    return value;
                }
                else{
                    double a1=0.0,a2=0.0,deno=0.0;
                    deno=t[i+k]-t[i];
                    if(deno==0){
                        a1=0;
                    }
                    else{
                        a1=(x-t[i])/deno;
                    }
                    deno=t[i+k+1]-t[i+1];
                    if(deno==0){
                        a2=0;
                    }
                    else{
                        a2=(t[i+k+1]-x)/deno;
                        
                    }
                    value1=a1*getBspline_value(x, i, k-1, t);
                    value2=a2*getBspline_value(x, i+1, k-1, t);
                    value=value1+value2;
                }
                return value;
            }
            else{
                cout<<"ERROR"<<endl;
                return 0;
            }
        }
    }*/
};


#endif /* ppform_h */
